#include <fstream>
#include <iostream>
#include <iomanip>
#include "TFile.h"
#include "TTree.h"
#include "TLine.h"
#include "TSystem.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TArc.h"
#include "TEllipse.h"
#include "TClonesArray.h"

using namespace std;

void neutron(){

  int NO = 890;

  int ret = gSystem->Load( "../Objects/libObjectClass.so" );
  if( ret!=0 ){
    cout << "Failed to Load Library: ret = " << ret << endl;
    return;
  }

  TFile *fread=new TFile("/scratch/shinichi/OutputFiles/AmBe.root");
  //  TFile *fread=new TFile("/scratch/shinichi/OutputFiles/9000neutron.root");
  if(!fread){ cout<<"No Such File"<<endl; return;}

  TTree *tree=(TTree*)fread->Get("tree");
  
  TClonesArray *JPMTa=new TClonesArray("PMT");
  TClonesArray *MPMT1a=new TClonesArray("PMT");
  TClonesArray *MPMT2a=new TClonesArray("PMT");
  TClonesArray *MPMT3a=new TClonesArray("PMT");
  TClonesArray *MPMT4a=new TClonesArray("PMT");
  TClonesArray *Particlea=new TClonesArray("Particle");
  TClonesArray *Neutrona=new TClonesArray("Particle");
  TClonesArray *Argona=new TClonesArray("Particle");
  TClonesArray *Gammaa=new TClonesArray("Particle");
  int nPhoton;
  int nElectron;

  tree->GetBranch("JPMTb")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT1b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT2b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT3b")->SetAutoDelete(kFALSE);
  tree->GetBranch("MPMT4b")->SetAutoDelete(kFALSE);
  tree->GetBranch("Particleb")->SetAutoDelete(kFALSE);
  tree->SetBranchAddress("JPMTb",&JPMTa);
  tree->SetBranchAddress("MPMT1b",&MPMT1a);
  tree->SetBranchAddress("MPMT2b",&MPMT2a);
  tree->SetBranchAddress("MPMT3b",&MPMT3a);
  tree->SetBranchAddress("MPMT4b",&MPMT4a);
  tree->SetBranchAddress("Particleb",&Particlea);
  tree->SetBranchAddress("Neutronb",&Neutrona);
  tree->SetBranchAddress("Argonb",&Argona);
  tree->SetBranchAddress("Gammab",&Gammaa);
  tree->SetBranchAddress("nPhoton",&nPhoton);
  tree->SetBranchAddress("nElectron",&nElectron);

  int nev = tree->GetEntries();
  cout<<"NEV = "<<nev<<endl;

  TGraph *NeutronXY[32];
  TGraph *NeutronXZ[32];
  TGraph *ArgonXY[32];
  TGraph *ArgonXZ[32];
  TGraph *GammaXY[32];
  TGraph *GammaXZ[32];
  for(int i=0;i<32;i++){
    NeutronXY[i]=new TGraph();;
    NeutronXZ[i]=new TGraph();;
    ArgonXY[i]=new TGraph();;
    ArgonXZ[i]=new TGraph();;
    GammaXY[i]=new TGraph();
    GammaXZ[i]=new TGraph();
  }


  TH1F *hNPhoton[8];
  TH1F *hNElectron[8];
  for(int i=0;i<8;i++){
    hNPhoton[i]=new TH1F(Form("hPhoton%d",i),Form("hPhoton%d",i),1000,0,100000);
    hNElectron[i]=new TH1F(Form("hElectron%d",i),Form("hElectron%d",i),500,0,50000);
  }

  int NNEUTRON,NARGON,NGAMMA;
  int Counter=0;
  int cESIS=0;
  int cESCP=0;
  int cISCP=0;
  int cES=0;
  int cIS=0;
  int cCP=0;
  int cTRANS=0;

  for(int iev=0;iev<nev;iev++){
    JPMTa->Clear();
    MPMT1a->Clear();
    MPMT2a->Clear();
    MPMT3a->Clear();
    MPMT4a->Clear();
    Particlea->Clear();
    Neutrona->Clear();
    Argona->Clear();
    Gammaa->Clear();

    tree->GetEntry(iev);

    bool IsElastic    = false;
    bool IsInElastic  = false;
    bool IsCapture    = false;

    //----- Neutron
    int ineutron = 0;
    int trackid = -1;
    int istep = 0;
    for(int i=0;i<Neutrona->GetEntriesFast();i++){
      Particle *neutron = (Particle*)Neutrona->At(i);
      if(trackid!=neutron->TrackID){
	trackid = neutron->TrackID;
	istep = 0;
	ineutron++;
	if( iev==NO ){
	  NeutronXY[ineutron-1]->SetPoint(istep,neutron->xpre,neutron->ypre);
	  NeutronXZ[ineutron-1]->SetPoint(istep,neutron->xpre,neutron->zpre);
	}
      }
      if( iev==NO ){
	//	cout<<istep<<" kine = "<<neutron->ekin<<endl;
	//cout<<istep<<neutron->processname<<endl;
	NeutronXY[ineutron-1]->SetPoint(istep+1,neutron->xpos,neutron->ypos);
	NeutronXZ[ineutron-1]->SetPoint(istep+1,neutron->xpos,neutron->zpos);
      }

      if     ( neutron->processname=="HadronElastic" )    IsElastic=true;
      else if( neutron->processname=="inelastic"     )  IsInElastic=true;
      else if( neutron->processname=="nCapture"      )    IsCapture=true;
      else if( neutron->processname!="Transportation") cout<<"ProcName = "<<neutron->processname<<endl;
      istep++;
    }
    if( iev==NO ){
      NNEUTRON = ineutron;
      cout<<"nNeutron = "<<NNEUTRON<<endl;
    }

    //----- Argon
    int iargon = 0;
    trackid = -1;
    istep = 0;
    for(int i=0;i<Argona->GetEntriesFast();i++){
      Particle *argon = (Particle*)Argona->At(i);
      if(trackid!=argon->TrackID){
	trackid = argon->TrackID;
	istep = 0;
	iargon++;
	if( iev==NO ){
	  ArgonXY[iargon-1]->SetPoint(istep,argon->xpre,argon->ypre);
	  ArgonXZ[iargon-1]->SetPoint(istep,argon->xpre,argon->zpre);
	}
      }
      if( iev==NO ){
	ArgonXY[iargon-1]->SetPoint(istep+1,argon->xpos,argon->ypos);
	ArgonXZ[iargon-1]->SetPoint(istep+1,argon->xpos,argon->zpos);
      }
      istep++;
    }
    if( iev==NO ){
      NARGON = iargon;
      cout<<"nArgon = "<<NARGON<<endl;
    }

    //----- Gamma
    int igamma = 0;
    trackid=-1;
    istep = 0;
    for(int i=0;i<Gammaa->GetEntriesFast();i++){
      Particle *gamma = (Particle*)Gammaa->At(i);
      if(trackid!=gamma->TrackID){
	trackid = gamma->TrackID;
	istep = 0;
	igamma++;
	if( iev==NO ){
	  GammaXY[igamma-1]->SetPoint(istep,gamma->xpre,gamma->ypre);
	  GammaXZ[igamma-1]->SetPoint(istep,gamma->xpre,gamma->zpre);
	}
      }
      if( iev==NO ){
	//	cout<<istep<<" edep = "<<gamma->edep<<endl;
	GammaXY[igamma-1]->SetPoint(istep+1,gamma->xpos,gamma->ypos);
	GammaXZ[igamma-1]->SetPoint(istep+1,gamma->xpos,gamma->zpos);
      }
      istep++;
    }
    if( iev==NO ){
      NGAMMA = igamma;
      cout<<"nGamma = "<<NGAMMA<<endl;
      cout<<"nPhoton = "<<nPhoton<<endl;
    }

    //-----
    hNPhoton[0]->Fill(nPhoton);
    hNElectron[0]->Fill(nElectron);
    Counter++;
    if( IsElastic&&IsInElastic ){
      hNPhoton[2]->Fill(nPhoton);
      hNElectron[2]->Fill(nElectron);
      cESIS++;
    }
    else if( IsElastic&&IsCapture ){
      hNPhoton[6]->Fill(nPhoton);
      hNElectron[6]->Fill(nElectron);
      cESCP++;
    }
    else if( IsInElastic&&IsCapture ){
      hNPhoton[7]->Fill(nPhoton);
      hNElectron[7]->Fill(nElectron);
      cISCP++;
    }
    else if( IsElastic ){
      hNPhoton[1]->Fill(nPhoton);
      hNElectron[1]->Fill(nElectron);
      cES++;
      if(nPhoton>0) cout<<"iev = "<<iev<<"/nphoton = "<<nPhoton<<endl;
      if(nElectron>0) cout<<"iev = "<<iev<<"/nelectron = "<<nElectron<<endl;
    }
    else if( IsInElastic ){
      hNPhoton[4]->Fill(nPhoton);
      hNElectron[4]->Fill(nElectron);
      cIS++;
    }
    else if( IsCapture ){
      hNPhoton[5]->Fill(nPhoton);
      hNElectron[5]->Fill(nElectron);
      cCP++;
    }
    else{
      hNPhoton[3]->Fill(nPhoton);
      hNElectron[3]->Fill(nElectron);
      cTRANS++;
    }
    
  }

  TCanvas *c1=new TCanvas("c1","c1",1200,600);
  c1->SetGrid();
  c1->Divide(2,1);
  c1->cd(1);
  TH1 *frame1=gPad->DrawFrame(-180,-180,180,180);
  frame1->SetXTitle("X (mm)");
  frame1->SetYTitle("Y (mm)");
  frame1->Draw();
  c1->cd(2);
  TH1 *frame2=gPad->DrawFrame(-180,-650,180,-300);
  frame2->SetXTitle("X (mm)");
  frame2->SetYTitle("Z (mm)");
  frame2->Draw();
  TLegend *leg1=new TLegend(0.6,0.2,0.85,0.35);
  leg1->SetFillStyle(0);

  for(int igamma=0;igamma<NGAMMA;igamma++){
    c1->cd(1);
    GammaXY[igamma]->SetLineColor(4);
    GammaXY[igamma]->SetMarkerColor(4);
    GammaXY[igamma]->SetMarkerStyle(21);
    GammaXY[igamma]->SetMarkerSize(1);
    GammaXY[igamma]->Draw("pl");
    c1->cd(2);
    GammaXZ[igamma]->SetLineColor(4);
    GammaXZ[igamma]->SetMarkerColor(4);
    GammaXZ[igamma]->SetMarkerStyle(21);
    GammaXZ[igamma]->SetMarkerSize(1);
    GammaXZ[igamma]->Draw("pl");
  }
  TLegendEntry *le=leg1->AddEntry(GammaXY[0],"Gamma","lp");
  le->SetTextColor(4);
  //
  for(int ineutron=0;ineutron<NNEUTRON;ineutron++){
    c1->cd(1);
    NeutronXY[ineutron]->SetLineColor(2);
    NeutronXY[ineutron]->SetMarkerColor(2);
    NeutronXY[ineutron]->SetMarkerStyle(21);
    NeutronXY[ineutron]->SetMarkerSize(1);
    NeutronXY[ineutron]->Draw("pl");
    c1->cd(2);
    NeutronXZ[ineutron]->SetLineColor(2);
    NeutronXZ[ineutron]->SetMarkerColor(2);
    NeutronXZ[ineutron]->SetMarkerStyle(21);
    NeutronXZ[ineutron]->SetMarkerSize(1);
    NeutronXZ[ineutron]->Draw("pl");
  }
  TLegendEntry *le=leg1->AddEntry(NeutronXY[0],"Neutron","lp");
  le->SetTextColor(2);
  //
  for(int iargon=0;iargon<NARGON;iargon++){
    c1->cd(1);
    ArgonXY[iargon]->SetLineColor(3);
    ArgonXY[iargon]->SetMarkerColor(3);
    ArgonXY[iargon]->SetMarkerStyle(21);
    ArgonXY[iargon]->SetMarkerSize(1);
    ArgonXY[iargon]->Draw("pl");
    c1->cd(2);
    ArgonXZ[iargon]->SetLineColor(3);
    ArgonXZ[iargon]->SetMarkerColor(3);
    ArgonXZ[iargon]->SetMarkerStyle(21);
    ArgonXZ[iargon]->SetMarkerSize(1);
    ArgonXZ[iargon]->Draw("pl");
  }
  TLegendEntry *le=leg1->AddEntry(ArgonXY[0],"Argon","lp");
  le->SetTextColor(3);
  //
  TArc *aInner=new TArc(0.,0.,133.,45.,315.);
  TLine *lInner=new TLine(133./sqrt(2),-133./sqrt(2),133./sqrt(2),133./sqrt(2));
  aInner->SetFillStyle(0);
  aInner->SetNoEdges(1);
  aInner->Draw();
  lInner->Draw();
  TLine *lFid1=new TLine(-133,-360,-133,-550);
  TLine *lFid2=new TLine(-133,-550,  94,-550);
  TLine *lFid3=new TLine(  94,-550,  94,-360);
  TLine *lFid4=new TLine(  94,-360,-133,-360);

  c1->cd(1);
  aInner->Draw();
  lInner->Draw();
  leg1->Draw();
  c1->cd(2);
  lFid1->Draw();
  lFid2->Draw();
  lFid3->Draw();
  lFid4->Draw();
  leg1->Draw();

  TCanvas *c2=new TCanvas("c2","c2",1200,600);
  c2->Divide(2,1);
  c2->cd(1);
  c2_1->SetGrid();
  //  c2_1->SetLogx();
  c2_1->SetLogy();

  TH1 *frame2_1=gPad->DrawFrame(0,0.1,10000,1e4);
  frame2_1->SetXTitle("# of Photon");
  frame2_1->SetYTitle("NEV");
  frame2_1->Draw();
  int color[8]={1,2,3,4,6,7,95,51};
  char name[8][32]={"All","Elastic","ES + IS","Transport","Inelastic","Capture","ES + CP","IS + CP"};

  TLegend *leg2_1=new TLegend(0.6,0.6,0.9,0.9);
  for(int i=0;i<8;i++){
    hNPhoton[i]->SetLineColor(color[i]);
    hNPhoton[i]->Draw("same");
    TLegendEntry *le=leg2_1->AddEntry(hNPhoton[i],name[i],"l");
    le->SetTextColor(color[i]);
  }
  leg2_1->SetFillStyle(0);
  leg2_1->Draw();

  c2->cd(2);
  c2_2->SetGrid();
  //  c2_2->SetLogx();
  c2_2->SetLogy();
  TH1 *frame2_2=gPad->DrawFrame(0,0.1,5000,1e4);
  frame2_2->SetXTitle("# of Electron");
  frame2_2->SetYTitle("NEV");
  frame2_2->Draw();
  TLegend *leg2_2=new TLegend(0.6,0.6,0.9,0.9);
  for(int i=0;i<8;i++){
    hNElectron[i]->SetLineColor(color[i]);
    hNElectron[i]->Draw("same");
    TLegendEntry *le=leg2_2->AddEntry(hNElectron[i],name[i],"l");
    le->SetTextColor(color[i]);
  }
  leg2_2->SetFillStyle(0);
  leg2_2->Draw();


  /*
  c2->cd(2);
  c2_2->SetGrid();
  c2_2->SetLogx();
  c2_2->SetLogy();
  TH1 *frame2_2=gPad->DrawFrame(10000,0.1,100000,1e4);
  frame2_2->Draw();
  int color[5]={1,2,3,4,6};
  for(int i=0;i<5;i++){
    hNPhoton[4-i]->SetLineColor(color[4-i]);
    if(i==1)
    hNPhoton[4-i]->Draw("same");
  }
  */

  cout<<"  All = "<<Counter<<endl;
  cout<<"   ES = "<<cES<<endl;
  cout<<"   IS = "<<cIS<<endl;
  cout<<"ES+IS = "<<cESIS<<endl;
  cout<<"ES+CP = "<<cESCP<<endl;
  cout<<"IS+CP = "<<cISCP<<endl;
  cout<<"   CP = "<<cCP<<endl;
  cout<<"Trans = "<<cTRANS<<endl;



}
